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We present an efficient Monte Carlo algorithm for determining the density of states which is based 
on the statistics of transition probabilities between states. By measuring the infinite temperature 
transition probabilities — that is, the probabilities associated with move proposal only — we are able 
to extract excellent estimates of the density of states. When this estimator is used in conjunction 
with a Wang-Landau sampling scheme [F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 
(2001)], we quickly achieve uniform sampling of macrostates (e.g., energies) and systematically refine 
the calculated density of states. This approach requires only potential energy evaluations, continues 
to improve the statistical quality of its results as the simulation time is extended, and is applicable 
to both lattice and continuum systems. We test the algorithm on the Lennard- Jones liquid and 
demonstrate good statistical convergence properties. 



I. INTRODUCTION 

In the fifty years since it was first introduced 0, the 
Monte Carlo (MC) method has become one of the pri- 
mary tools for the simulation of materials and the pre- 
diction of their properties. MC simulations are widely 
used and a significant number of methodological im- 
provements now exist for efficient and accurate studies of 
systems which are otherwise intractable with the tradi- 
tional Metropolis sampling |2|.|3| Examples include high- 
density or low-temperature liquids, network-forming flu- 
ids, polymers, and proteins, for which the simulated sys- 
tem is prone to becoming trapped in potential energy 
minima for large numbers of simulation steps. 

A number of these extended MC methods solve the 
ergodicity problem by forcing a broad sampling of phase 
space; included in this category are parallel tempering 
multicanonical methods J5J, and the more recent Wang- 
Landau (WL) algorithm |fg. The latter two achieve broad 
phase space sampling by attempting to produce a uni- 
form distribution of one or more macroscopic observables, 
such as number of particles, volume, and potential en- 
ergy. To arrive at such a distribution, each microstate 
of the system must be sampled with a probability in- 
versely proportional to the degeneracy of the correspond- 
ing macroscopic observable. Thus whether implicitly (as 
part of a state weight) or explicitly (through direct cal- 
culation), these methods must determine the density of 
states, f2. Of course, fi contains a complete thermody- 
namic description of the system via Boltzmann's equa- 
tion, 

S(N, V,E) = klnQ(N,V,E). (1.1) 

Clearly, the density of states is not known prior to the 
simulation but must be successively approximated dur- 
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ing its course. The multicanonical approach iteratively 
updates its estimate for weights (which include f2) from 
histogram results of several pre-production runs; once 
the density of states is reasonably converged, it and his- 
togram data from a final long simulation interval are 
"reweighted" to determine properties at the state condi- 
tions of interest In the WL scheme, microstates are 
explicitly sampled with the probability fl^ 1 and the den- 
sity of states is modified at every simulation step [|| 0] ; 
for each visited state, the corresponding value of £1 is 
scaled by the modification factor, /. In this way, the 
the density of states is designed to converge to its cor- 
rect value when macroscopic properties are sampled uni- 
formly. However, to eliminate / fluctuations in the final 
result, a schedule must be imposed in which / changes 
from an initially large value (say, / ~ e) to nearly unity. 
Often, this schedule is implemented by maintaining a his- 
togram of visited states, and when the bin with the fewest 
entries is no less than 80% of the average bin (i.e., when 
the histogram is sufficiently flat), / is changed according 
to /now = \/ /old and the histogram is reset. 

The WL algorithm has received much attention lately, 
perhaps due to its straightforward implementation and 
wide applicability, and a number of methodological vari- 
ants have been proposed, including multibondic |8j, N- 
fold way 0, and quantum 01 extensions. Of relevance 
to continuum systems, Yan and de Pablo noticed two 
deficiencies of the WL method: estimates of the density 
of states reach a limiting statistical accuracy which is 
not improved with further MC steps, and the large num- 
ber of configurations generated towards the end of the 
simulation make only a small contribution to the calcu- 
lated density of states. They proposed two variants in 
which estimates of temperature improve the convergence 
of the density of states relative to the original method 
[Tl| . Both of these, however, must be formulated specifi- 
cally for the relevant macroscopic parameter (the authors 
originally discussed only energy in this regard). Further- 
more, their configurational temperature method requires 
first and second derivatives of the potential energy, which 
limits it to systems with continuous potentials and, prac- 
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tically speaking, excludes those for which such calcula- 
tions are computationally expensive. 

Here we present a MC implementation for directly de- 
termining the density of states of a system, which we be- 
lieve is superior to original WL algorithm while retaining 
much of its generality. A reformulation of the Wang- 
Landau approach, this implementation is based on the 
ideas of transition matrix MC methods: the calculations 
utilize a record of proposed transitions between, rather 
than visits to, various macrostates 0, 0, 0, 0, 0, 0, 
H S El El El E| . Our method builds on the ideas 
of Ref. 11 in that we supplement the normal WL sam- 
pling scheme with this additional data collection, which 
ultimately produces density of states estimates of higher 
statistical quality than the WL estimates. In this paper 
we show that our method is: (1) general — applicable to 
lattice and continuum systems; (2) robust — continuously 
improves its estimate for the density of states, regard- 
less of a particular modification factor schedule; (3) sim- 
ple — requires only potential energy evaluations; (4) ef- 
ficient — uses all collected information during the course 
of a simulation; and (5) accurate — produces density of 
states estimates of higher statistical quality than the WL 
method for a given number of simulation steps. 

In what follows we first outline the theoretical ba- 
sis for transition matrix calculations and review a pre- 
vious method based wholly on transition probabilities 
[TH EH EH E3- We then present our implementa- 
tion, which combines transition matrix ideas with Wang- 
Landau sampling, and discuss its advantages relative to 
both the WL and transition matrix formulations alone. 
Finally, we proceed to discuss results for a set of test 
cases. 



II. THEORETICAL BASIS 

We begin by considering the transition probability 
of moving from macrostate I to J given that the sys- 
tem is currently in /, which we denote by T(I — > J). 
This derivation applies generally to any well-defined 
macrostate; in this work, however, we use / and J to de- 
note energy levels. (We refer the reader to Refs. El and 
l24l which demonstrate the use of transition probabilities 
with a different macrostate variable.) The expression for 
T(I — > J) in terms of the corresponding microstate tran- 
sition probabilities is then 

T(I - J) = n(I)- 1 £ J2 t(i -> 3) Sa 5, j (2.1) 

* i 

where the sums extend over all microstates, f2(7) gives 
the degeneracy of macrostate I, t(i — > j) is the probabil- 
ity that the system will make a move from microstate i 
to j given that it is currently in state i, and Su (equiva- 
lently, 6jj) is defined such that 

6,1 = { otherwise. ^ 



That is, Sn is one if and only if microstate i is a member 
of macrostate I. Implicit in Eg. 12. II is the assumption of 
equal a priori probability of microstates; the macrostate 
transition probability is simply the microcanonical av- 
erage of all corresponding microstate transition proba- 
bilities. Accordingly, the following results apply to any 
simulation procedure which generates microcanonical av- 
erages according to Eq. 12.11 By dividing Eq. 12.11 by 
the corresponding expression for the reverse transition, 
we obtain an important relation between the transition 
probabilities and density of states: 

t(i —»«/) _ n( J) ( E t E 3 t(i -> 3) Su Sjj \ 

nj^i) m [ Ej tu - ») ^ y • 1 ■ ' 

The microstate transition probability t(i — > j) appear- 
ing in these expressions is the product of two terms: 
a(i — > j), the probability of proposing a transition from 
i to j, dependent only on the type of MC move em- 
ployed; and Pncc(i — > j), the probability of accepting a 
proposed move from i to j, dependent on the state con- 
ditions of the particular statistical ensemble. We con- 
sider the special case in which moves are strictly ac- 
cepted, P acc (i ~^ j) — 1) which corresponds to canonical 
Metropolis-style sampling at infinite temperature |2fij . 
Consequently, we denote (I — > J) as that correspond- 
ing to strict move acceptance and will refer to it as an 
infinite temperature transition probability. We note that 
-Pacc(* ~~ * j) = 1 is not necessarily the infinite-T limit 
for all ensembles and classes of simulation moves; how- 
ever, in identifying T 00 (I — > J) with the strict acceptance 
condition, this derivation is in fact quite general, regard- 
less of the particular move type or ensemble used. The 
relationship among the infinite-temperature macrostate 
transition probabilities is then 

r^ji —>«/) _ njj) ( E, E 3 oft ^ i) Su SjA 

^ ( J - I) m { Ei E, <3 - i) Su hi)' [ ' 1 

Given the symmetry properties of a particular move type 
[i.e., the ratio a(i — > j)/ct(J — » i)], this equation permits 
an estimate of the density of states from knowledge of 
macrostate transition probabilities. In the case of sym- 
metric moves with a(i —* j) = a(J — > i), the term in 
parenthesis in Eq. 12.41 drops out completely; examples 
of such moves include single-spin flips and single-particle 
displacements. It is important to consider T 00 (J — > J) 
as the probability that a move will be proposed from / 
to J. Thus, the collection of Too (I — > J) quantities gives 
a measure of the connectivity of macrostates for a given 
move type. Previous studies have demonstrated that this 
type of information is statistically superior to histogram 
data 0, 0, E3 ) an d its use is central to our method. 

Though the preceding derivation entails discrete 
macrostates, the extension to off-lattice systems is con- 
ceptually identical and thus relatively straightforward 
|22j . Practically speaking, however, one must impose 
a discretization scheme on any continuous macroscopic 
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variables for the measurement of transition probabilities 
in a simulation. Because of this discretization, the mi- 
crostate sampling scheme must also be discretized. For 
example, each configuration's energy must be assigned an 
effective value, the average energy of the corresponding 
bin, for use in the acceptance criterion. Lack of congru- 
ence with the transition matrix discretization violates the 
microcanonical average in Eq. 12.11 and results in system- 
atic errors in the density of states. 



III. ALGORITHMIC DETAILS 

To calculate Q for a given system, two algorithmic 
ingredients are necessary, as discussed by Oliveira, et. 
al. (i) a good way to measure macrostate degen- 

eracies, termed an "estimator," and (ii) broad sampling 
of phase space such that all macrostates of interest are 
visited. The first condition is met by measurement of 
macrostate transition probabilities and use of Eq. 12.41 
One can estimate infinite-T transition probabilities by 
recording move proposal statistics during the course of a 
simulation, even though the acceptance of the moves may 
not correspond with the infinite-T condition [2(1 l2ll l22j | . 
In other words, the T 00 (/ — > J) values can be determined 
using any arbitrary sampling scheme for P acc {i —* j) since 
they depend only on the proposal probabilities contained 
in a(i — ► j). The procedure involves adding 1 to a matrix 
entry C(I, J) every time a move from macrostate I to J 
is proposed. The estimate for the transition probability 
is 



T 00 (I^J) = C(I,J)/J2C(I,K) 



(3.1) 



K 



where sum extends over all macrostates and the tilde 
indicates the estimate. Once the infinite-T transition 
probabilities are known, the density of states can be cal- 
culated from Eq. 12.41 In general, this is an over-specified 
problem since in a system of N macrostates with N un- 
known values of O (more accurately, N — 1 unknown 
relative values), there are N(N — l)/2 such equations. 
In solving for the density of states given the transition 
probabilities, the easiest approach considers only neigh- 
boring states for N — 1 instances of Eq. 12.41 However, 
one can consider all N(N — l)/2 equations by minimiz- 
ing the variance among the predicted values of lnfi, as 
described by Ref. H3 Denoting S(I) = \nfl(I) and 
H(I) = ^2 K C(I,K), one minimizes the total variance 
with respect to the values S(I), where the variance is 
defined as 



°tot = ^2 



[S(I) - S(J) + lnTooOT - JVTooG/ -> I)] 2 



i.j 



(T 



IJ 



(3.2) 



with 



ij 



Here, we use the number of entries in the C matrix to 
determine the weight of each estimate in the minimiza- 
tion procedure, assuming that var [C{I, J)] cx C(I, J). In 
the solution to these equations, one value of the density 
of states must be fixed and the remaining can be solved 
by matrix inversion. 

The second condition for calculation of the density of 
states, broad phase space sampling, can be accomplished 
by using a uniform ensemble in which all macrostates are 
equiprobable. In such an ensemble, the probability of a 
given microstate is inversely proportional to the density 
of states of the corresponding macrostate. Accordingly, 
the acceptance criterion for symmetric moves is 



■Pacc(« -> j) = min 



1 wy\ 

>n(j)J 



(3.4) 



Of course, the density of states which comes into play in 
this expression is not known a priori. A straightforward 
implementation of transition probabilities is possible by 
combining Eq. 12.41 with this acceptance criterion : 



-P aC c(« -> j) = min 1, 



Tqo(J^7) 



(3.5) 



C(I, J)- 1 + H(I)- 1 + C(J, I)- 1 + H(J)-\ (3.3) 



With Eas. 13.11 and !3. 51 one can immediately construct a 
complete MC procedure, an idea originally proposed by 
Oliveira, et. al. [Ta and later generalized by Wang and 
coworkers [2(], l2lL l22j|. The simulation is started with 
C(I, J) = for all / and J. Moves are then proposed 
and accepted/rejected based on the acceptance criterion 
in Eq. 13.51 and after each such step, the C matrix is up- 
dated to reflect the proposed move, regardless of whether 
or not it was accepted. In this fashion, the approxi- 
mations Too (I —> J) grow more accurate and the dis- 
tribution of sampled macrostates becomes increasingly 
more uniform. Here, it is important to note that even 
though the actual acceptance probability is not that in 
the infinite-T case, the estimates for T 00 (7 — > J) are con- 
structed as if P&c C (i — » j) — 1 due to the way the C 
matrix is updated. The density of states is never directly 
referenced in this scheme. Instead, one is concerned only 
with the infinite-T transition probabilities, and therefore 
need only maintain the C matrix. At the end of the 
simulation or at any intermediate point, an estimate for 
density of states can be determined using Eq. 12.41 This 
estimate improves as the number of steps in the simula- 
tion increases, and furthermore, the quality of the esti- 
mate can be assessed by the number of entries in each 
C(I,J) 01 . In this paper, we term this approach the 
transition matrix (TM) method. 

The procedure just described is essentially an iterative 
scheme for the solution of a set of nonlinear equations 
describing the density of states |2£J, l2ll |22j • In general, 
this is a problem for which convergence is not immedi- 
ately obvious nor guaranteed. In particular, problems 
arise with the acceptance criterion in Eq. 13.51 early on 
in the simulation when numerous zeros exist in the C 
matrix; an estimate T 00 (/ — > J) is ill-defined if either 
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C(I, J) — or J2k K) — 0. Consequently, for runs 
which probe large ranges of O, the time required for all 
macrostates to be initially visited can become extraor- 
dinarily long. In contrast, the Wang-Landau approach 
forces the complete set of macrostates to be visited very 
quickly at the beginning of the simulation, but lacks the 
accuracy and continuous improvement in calculating the 
density of states that the TM method exhibits. 

Fortunately, one can combine the strengths of both 
algorithms, an implementation that we denote as the 
WL-TM method. In this approach, we use the origi- 
nal acceptance criterion in Eq. 13.41 and update the den- 
sity of states every move as in the WL method, using 
£1(1) *— / x £1(1) for the state at the end of each MC 
step. In addition, we start with a zeroed C matrix, as in 
the TM method, and record all proposed transitions dur- 
ing the entire schedule of modification factors (i.e., the C 
matrix is never re-zeroed, and therefore contains informa- 
tion from the complete duration of the simulation) . Peri- 
odically, a completely new density of states is generated 
from the infinite temperature transition probabilities via 
Eqs. 12.41 and 13.11 — we call this "refreshing" the density 
of states. In this way, the WL scheme quickly enforces a 
broad distribution of macrostates, but the TM refreshing 
accelerates the convergence of the density of states and 
hence speeds the modification factor schedule. In the lat- 
ter portions of a simulation, the modification factor is so 
close to one that its effects on £l are negligible. However, 
the refreshing continuously improves the calculation since 
the statistics of the C matrix grow better with the num- 
ber of simulation steps. At this point, it is even feasible 
to switch to the bare TM scheme described by the ac- 
ceptance criterion in Eq. 13.51 Therefore, the distinction 
between our approach and the TM method is that we 
use the WL sampling scheme to guarantee convergence of 
the infinite-T transition probabilities, whereas the latter 
method has no such guarantee and might requir e a good 
initial estimate for its transition probabilities 20, 2lll22|. 

It is important to note that the modification factor 
schedule in this hybrid scheme might be modified from 
the original WL approach. For the Lennard-Jones sys- 
tem, we have found that the 80% "flat" histogram re- 
quirement for changes to / can be relaxed significantly to 
the mere requirement that each histogram bin have only 
one entry. Furthermore, updates to the modification fac- 
tor can be more dramatic, such as In / <— 0.1 In/. These 
changes minimize the time that detailed balance is not 
obeyed, during which entries to the C matrix might be 
biased due to violation of the microcanonical average in 
Eq. 12.11 As In / approaches zero, updates to the infinite 
temperature transition probabilities grow more accurate. 
That is, the faster In / approaches a near-zero value, the 
more accurately the density of states estimate will be 
calculated from the transition probabilities. One might 
consider postponing entries in the C matrix and the re- 
freshing of the density of states until the first few stages 
of the modification factor schedule have passed, during 
which / 3> 1. It is also possible to refresh only once, at 



10 1 r 



c<3 

.£ 10 u t 
o 

10 k 

to 
0) 

< 10" 2 fe- 



10"- 



TTTl I 



I | 1 1 I I I 1 1 1 1 



♦ * * 

□ 



□ 





LJ (WL) 




♦ 


LJ (TM) 




□ 


LJ (WL-TM) 




• 


Ising (TM) 





eft 



9 



m 



9 9 



m -. 



ul_ 



10 10 10 

Millions MC steps 



10 J 



FIG. 1: Error analysis of the dimensionless entropy, S = InQ, 
as calculated by the Wang-Landau (WL), transition-matrix 
(TM), and combination (WL-TM) methods for the 110 par- 
ticle Lennard-Jones (LJ) and 32x32 Ising systems. Details of 
each method are described in the text. 



the end of the simulation, which is equivalent to the orig- 
inal WL algorithm with additional data collection (the 
transition matrix) . These options impart some flexibility 
to the WL-TM approach, to be tailored to the specific 
system of interest (for the Lennard-Jones system, which 
converges relatively fast, we have found little difference 
among them in our tests). 



IV. TEST CASES 

We have tested these algorithms on the 110-particle 
Lennard-Jones system at reduced density p — 0.88. We 
cut the potential at 2.5a and apply the long-range tail 
correction and we tabulate the transition probabil- 
ities for the energy range — 700e to — 500e, discretized 
into 100 bins. We reject moves which would take the 
system out of this energy range and accordingly update 
the transition matrix by C(J, I) «— C(I, I) + 1 where I 
is the initial state; this has the same effect on the tran- 
sition probabilities and subsequent density of states cal- 
culation as if the range were not truncated. As in the 
usual implementation of these methods, we calculate the 
dimensionless entropy, S = \n£l, rather than the density 
of states itself. 

We perform six independent runs for each simulation 
and measure the standard deviation of corresponding en- 
tropy bins among them; the statistical error is then mea- 
sured as the average of this standard deviation over all 
the bins. (This approach fails to capture systematic er- 
rors; therefore, below we also discuss agreement with re- 
sults from traditional MC simulations.) For comparison, 
we also performed the TM method for the 32x32 Ising 
system, employing single-spin flips and tabulating tran- 
sition probabilities for all 1023 energy levels. The esti- 
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FIG. 2: Comparison of traditional methods (rilled symbols) 
with calculations from density of states generated with the 
combination (WL-TM) method and reweighted to the appro- 
priate state conditions (lines). 



FIG. 3: Errors for each method in the calculated isochoric 
heat capacity, plotted as a function of CPU time. For the 
TM method, results from the early portion of the simulation, 
during which the relevant energy levels have not yet been 
visited, are omitted. 



mate of statistical accuracy for the Ising system is the 
fractional error, 1 — |S'run/>S'exact|) averaged over each en- 
ergy level and three independent runs; the exact results 
are from the method of Ref . [27J • 

The results of these tests are shown in Fig.^ One can 
see that statistical errors in the WL method rapidly de- 
crease to a limiting value, beyond which additional MC 
steps do not improve the calculations, as also shown in 
Ref. 0. This is primarily due to the fact that changes 
in the modification factor eventually become too fast rel- 
ative to its effect on the density of states. For the TM 
method, errors are significant early on in the simulation 
due to the time needed to generate initial transition prob- 
ability estimates for each energy level. Once this range 
has been covered, the density of states is systematically 
refined with much greater precision than the WL method. 
(Our experience with the TM approach indicates that 
this initial lag time can grow significantly worse as the en- 
ergy range broadens, eventually making the TM method 
infeasible for large entropy gradients.) The WL-TM al- 
gorithm successfully overcomes this drawback and has 
good convergence properties in both the initial and final 
stages of the simulation. 

For the Lennard-Jones system, we have used a stan- 
dard reweighting procedure [2g to generate predictions 
for the potential energy from the density of states output 
of the WL-TM algorithm. In Fig. we show these re- 
sults alongside data taken from standard NVT MC sim- 
ulations [23. The results from the WL-TM algorithm, 
which are taken only 10 million steps into the simula- 
tion, are in excellent agreement with the results of the 
traditional method. In Fig. [31 we also show the statis- 
tical error in the calculated heat capacity at T = 1.0, 
which is measured among results from six independent 
runs. The heat capacity results of the WL-TM algo- 
rithm appear to have similar convergence behavior as the 



methods introduced by Yan and de Pablo 0] > indicating 
that for simple systems, it does not offer appreciable im- 
provement over these previous approaches. However, the 
current method achieves this level of performance with- 
out energy derivative information and while remaining 
readily extendable to any macrostate parameter, includ- 
ing arbitrarily defined ones (e.g., as in extended ensemble 
formulations) . 

In Fig. 01 we demonstrate the effects of system size and 
energy discretization on the WL-TM algorithm. Results 
are shown for the original system as well as for a 250 
particle system at the same density but with roughly the 
same energy discretization (~ 2e). Here, the statistical 
error in the dimensionless entropy is normalized by the 
number of particles, and therefore effectively measures 
a fractional error. Fig. 0] shows that the convergence of 
the two system sizes is comparable; the 250 particle case 
does not noticeably suffer from methodological perfor- 
mance losses. This emerges from the banded nature of 
the transition matrix: each energy level is effectively con- 
nected to only a few energy neighbors due to the small 
range of energy explored by single particle moves, which 
have been tuned for 50% acceptance. We also investigate 
the effects of energy discretization on the calculations by 
increasing the number of energy bins while keeping the 
system size constant at 110 particles. Fig. 0]reveals that 
fine-graining the energy range only weakly affects the sta- 
tistical errors, despite the fact that the number of entries 
in the C matrix (each for which an estimate must be 
generated) grows as the square of the number of bins. 
Presumably this occurs due to the banded nature of the 
matrix and the subsequent error minimization procedure, 
which takes into account all transition state pairs in de- 
termining the density of states. 

We have also investigated our approach on a more 
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FIG. 4: Effects of system size and energy discretization on the 
WL-TM method. All errors are measured for the dimension- 
less entropy, S, and are normalized by the system size. The 
upper points, which have been scaled by a factor of 100 for 
clarity, demonstrate the effects of increasing the energy dis- 
cretization for a fixed system size of 110 particles. The lower 
points correspond to system sizes of 110 and 250 particles at 
fixed density for a constant degree of discretization (~ 2e). 



challenging system, liquid silica as described by the Van 
Beest-Kramer-Van Stanten (BKS) potential [3(j. This 
system is one for which potential energy derivatives are 
expensive due to long-range interactions, which encour- 
ages methods employing only energy evaluations. The 
BKS potential treats silica as independent oxygen and 
silicon ions which interact pairwise through electrostatic 
and Buckingham-type repulsive/dispersive forces. Long- 
range Coulombic interactions are calculated with the 
Ewald method, using all /c-vectors with magnitude less 
than or equal to 5 x 2n/L and aL — 5.6. Short-range 
interactions are truncated and shifted at 5.5 A. We per- 
form the simulation for 100 oxygen and 50 silicon atoms 
at a density of 2200 kg/m 3 , and sample the energy range 
-5.5 to -5.2 MJ/mol, discretized into 200 bins. At this 
density, the dynamics of the silica system are known 
to experience a dramatic slowdown as the temperature 
is decreased, which traditionally requires long molecular 



dynamics trajectories to equilibrate [3l|. We run the pro- 
posed WL-TM algorithm for the small silica system for 
5 x 10 8 MC steps. For comparison, we conduct multiple 
NVE molecular dynamics (MD) simulations correspond- 
ing to 11 different temperatures, using velocity- Ver let in- 
tegration with a 1 fs timestep and propagating for 2 x 10 5 
timesteps. Fig. [5] reveals good agreement between the re- 
sults of the WL-TM algorithm for silica and those gener- 
ated from the MD simulations. This is merely a first step 
in demonstrating the usefulness of WL-type methods for 
the glassy silica system, and we leave a more detailed 
investigation for future work. 

V. CONCLUSION 

In summary, we have presented a method for calcula- 
tion of the density of states of a system which combines 
the good statistical accuracy of transition matrix estima- 
tors with the rapid broad sampling of phase space gener- 
ated by the Wang-Landau method. This implementation 
requires only potential energy evaluations, and therefore, 
is general to lattice and continuum systems. For the same 
reason and for the additional feature that the transition 
matrix preserves information from the complete simu- 
lation, the method is computationally efficient. Finally, 
from our test simulations, our implementation appears to 
have good convergence properties and allows rapid calcu- 
lation of the density of states as compared to the original 
Wang-Landau approach. The method appears encour- 
aging for studies of complex systems such as dense, su- 
percooled, or glassy liquids, for which gaining accurate 
results from reasonable simulation times has been chal- 
lenging. 
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